
#################################################
####### Replication code for "The Connections of Party Brokers," Journal of Politics
####### This file includes replication code for all analyses in the MAIN TEXT
###################################################

rm(list=ls())
library(foreign)
library(MASS)
library(tidyverse)
library(psycho)
library(stargazer)
library(margins)
library(dotwhisker) 
library(broom) 
library(modmarg)
library(xtable) 
library(stargazer)
library(apsrtable)




###################################
###### CALL IN DATA
###################################

load("BrierleyNathan_brokers.Rdata")

#Contains the following 4 objects:
#d: individual-level survey data. row = aspiring branch leader in 2018
#bl: d collapsed to the Branch Level ("bl"). row = polling station. variables such as connections up/down are the average across the incumbent branch leaders
#ps: d collapsed to the Position Level ("ps"). row = all contestants, for example, for secretary within branch X. variables such as connections up/down are for the incumbent for that position
#elections: is d subset to participants in actual contested elections for positions

#d.incrun: d subset only to incumbent branch executives who are recontesting (for use in subsequent analyses below)
d.incrun <- dplyr::filter(d, old_branch_position==1 & run_same_pos==1)
dim(d.incrun)


######################################################################
######################################################################
### CODE FOLLOWS ORDER OF PRESENTATION IN MAIN TEXT
######################################################################
######################################################################


###################################
###### Figure 3: Individual-level correlations between connections up and connections down 
###################################

#FIGURE 3
par(mar=c(5,5,3,3))
plot(jitter(d$cxn_down_percentage_correct_full), jitter(d$cxn_up_percentage_correct_full), axes=F, ylab="Connections up (proportion correct)", xlab="Connections down (proportion correct)", pch=19, col="white", xlim=c(0,1), cex.lab=1, ylim=c(0,1), cex.axis=1.2) 
     #main="Individual-level correlation\n(overall sample)")
points(jitter(d$cxn_down_percentage_correct_full), jitter(d$cxn_up_percentage_correct_full),  col="darkgrey", pch=20)
axis(1, cex.axis=1.2)
axis(2, cex.axis=1.2)
#lines(lowess(d$cxn_down_percentage_correct_full, d$cxn_up_percentage_correct_full), col="red") # lowess line (x,y)
abline(lm(d$cxn_up_percentage_correct_full~d$cxn_down_percentage_correct_full), col="dodgerblue2", lwd=2.5, lty="solid") # regression line (y~x) 
#abline(v=mean(d$cxn_up_percentage_correct_full), lwd=1.2, col="firebrick", lty="dashed")
#abline(h=mean(d$cxn_down_percentage_correct_full), lwd=1.2, col="firebrick", lty="dashed")
text(0.85, 0.9, "r = -0.003", cex=1.2)

#summary statistics in the paper:
summary(d$cxn_down_percentage_correct_full)
#   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#0.00000 0.04167 0.16667 0.21447 0.33333 1.00000 
sd(d$cxn_down_percentage_correct_full)
#0.2219237
summary(d$cxn_up_percentage_correct_full)
# Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
# 0.0000  0.0800  0.2000  0.2152  0.3200  0.8400 
sd(d$cxn_up_percentage_correct_full)
#0.1525035


#This is true in the full sample and within either rural or urban constituencies. 
summary(lm(cxn_up_percentage_correct_full ~ cxn_down_percentage_correct_full + urban, data=d))
summary(lm(cxn_up_percentage_correct_full ~ cxn_down_percentage_correct_full, data=d))

###################################
###### Figure 4: Predictors of connections down (left) and connections up (right).  
###################################


#first: select variables to standardize 
d_z1 <- select(d, age, asset_index, years_comm)
d_z2<- standardize(d_z1)
#second: select other relevant vars
d_z3 <- select(d, cxn_down_percentage_correct_full,cxn_up_percentage_correct_full,  female, local_eth_minority, bio_educ_three, assoc_community,  petty_trader,  farmer ,  inf_business ,  civil_servant ,  chief_relative ,  pol_relative ,  anon.ps.name)
#third: join the two datasets
d_z4 <- cbind(d_z2, d_z3)


### Run the regs on the standardized dataset (z_d)
m1.c.fe.ps <- lm(cxn_down_percentage_correct_full ~ female +  age +local_eth_minority + as.factor(bio_educ_three) + asset_index + years_comm + assoc_community+ petty_trader + farmer + inf_business + civil_servant + chief_relative + pol_relative + as.factor(anon.ps.name), data=d_z4)
summary(m1.c.fe.ps)

m2.c.fe.ps <- lm(cxn_up_percentage_correct_full ~  female +  age +local_eth_minority + as.factor(bio_educ_three) + asset_index   + years_comm  + assoc_community + petty_trader + farmer + inf_business + civil_servant + chief_relative + pol_relative + as.factor(anon.ps.name), data=d_z4)
summary(m2.c.fe.ps)


par(mfrow=c(1,2))

par(mar=c(5,11,3,2))
plot(coef(m1.c.fe.ps)[2:15], c(14:1), xlim=c(-0.1, 0.1), col="white", xlab="Standardized Coefficient Estimate", yaxt="n", ylab="", main="Predictors of Connections Down\n(proportion correct)")
abline(v=0, col="firebrick", lty="dashed", lwd=0.9)
axis(2, a=c(14:1), labels=c("Female","Age", "Local eth. minority", "Secondary education", "Tertiary education", "Assets", "Years in community", "Association member", "Petty trader", "Farmer", "Small businessperson", "Civil servant", "Relative of chief", "Relative of pol."), las=1)
segments(confint(m1.c.fe.ps)[2:15,1], c(14:1), confint(m1.c.fe.ps)[2:15,2], c(14:1), col="grey60", lwd=2)
segments(confint(m1.c.fe.ps, level=0.9)[2:15,1], c(14:1), confint(m1.c.fe.ps, level=0.9)[2:15,2], c(14:1), col="black", lwd=2)
points(coef(m1.c.fe.ps)[2:15], c(14:1), pch=16, col="dodgerblue2", cex=1.4)

par(mar=c(5,11,3,2))
plot(coef(m2.c.fe.ps)[2:15], c(14:1), xlim=c(-0.1, 0.1), col="white", xlab="Standardized Coefficient Estimate", yaxt="n", ylab="", main="Predictors of Connections Up\n(proportion correct)")
abline(v=0, col="firebrick", lty="dashed", lwd=0.9)
axis(2, a=c(14:1), labels=c("Female","Age", "Local eth. minority", "Secondary education", "Tertiary education", "Assets", "Years in community", "Association member", "Petty trader", "Farmer", "Small businessperson", "Civil servant", "Relative of chief", "Relative of pol."), las=1)
segments(confint(m2.c.fe.ps)[2:15,1], c(14:1), confint(m2.c.fe.ps)[2:15,2], c(14:1), col="grey60", lwd=2)
segments(confint(m2.c.fe.ps, level=0.9)[2:15,1], c(14:1), confint(m2.c.fe.ps, level=0.9)[2:15,2], c(14:1), col="black", lwd=2)
points(coef(m2.c.fe.ps)[2:15], c(14:1), pch=16, col="dodgerblue2", cex=1.4)

#cadetblue1

###################################
###### Figure 5: Correlations between connections down and individual- and community-level variables. 
###################################


cor(d$cxn_down_percentage_correct_full, d$chief_relative)
#0.2613782
cor(d$cxn_down_percentage_correct_full, d$years_comm)
#0.2745849
cor(d$cxn_down_percentage_correct_full, d$lives_outside_ps)
#-0.1127439
cor(d$cxn_down_percentage_correct_full[is.na(d$TOTAL.VOTERS.x)==F], d$TOTAL.VOTERS.x[is.na(d$TOTAL.VOTERS.x)==F])
#-0.3362391
cor(d$cxn_down_percentage_correct_full, d$farming_households_2km_p)
#0.6415739
cor(d$cxn_down_percentage_correct_full[is.na(d$closest_ps_km.x)==F], d$closest_ps_km.x[is.na(d$closest_ps_km.x)==F])
#0.4319898
cor(d$cxn_down_percentage_correct_full, d$urban)
#-0.4876391
cor(d$cxn_down_percentage_correct_full, d$wealth_index_2km)
#-0.63
cor(d$cxn_down_percentage_correct_full, d$ethfrac_upperlevel_2km)
#-0.3151215


## PLOTTING THE CORRELATION COEFS
def_y <- seq(from = 1, to = 9, by = 1)
def_x <- seq(from = -1, to = 1, by = 0.25) 

## Need to change plot area
par(mfrow=c(1,1))
par(mai=c(1.02,2.5,0.82,0.42))

plot(def_x, def_y, col="white", axes=F, ylab="", xlab="Correlation coefficient", cex.lab=1.2)
points(-0.63, 1, pch=19,col="dodgerblue2", cex=1.6) 
points(-0.49, 2, pch=19,col="dodgerblue2", cex=1.6) 
points(-0.34, 3, pch=19,col="dodgerblue2", cex=1.6) 
points(-0.32, 4, pch=19,col="dodgerblue2", cex=1.6) 
points(-0.11, 5, pch=19,col="dodgerblue2", cex=1.6) 
points(0.26, 6, pch=19,col="dodgerblue2", cex=1.6) 
points(0.27, 7, pch=19,col="dodgerblue2", cex=1.6) 
points(0.43, 8, pch=19,col="dodgerblue2", cex=1.6) 
points(0.64, 9, pch=19,col="dodgerblue2", cex=1.6) 
abline(v=0, lty=2, lwd=1.6)
axis(2, at=(1:9), labels=c(
  "Wealth index (PS)", 
  "Urban (PS)", 
  "# of registered voters (PS)", 
  "Ethnic frac. (PS)", 
  "Lives outside PS (IND)", 
  "Relative of chief (IND)", 
  "Years in community (IND)",
  "Distance to next PS (PS)",
  "% farming households (PS)"), las=2, cex.axis=1) 
axis(1, cex.axis=1.2)




###################################
###### Table 1: Summary of activities that branch leaders perform 
###################################

d_brok <- select(d, cpgn_brok_index, campaign_canvass,
                        campaign_orgrallies,
                        campaign_orgevents, 
                        campaign_gifts, 
                        campaign_orgtransport,
                        campaign_assistance, 
                        campaign_workwchief,
                        campaign_drivevoters,
                        campaign_jobs,
                        broker_up_down_index, 
                        helps_comm_w_party,
                        helps_comm_w_da,
                        helps_party_idcomm)
stargazer(d_brok, iqr=F)



###################################
###### Table 2: Mean connections down of aspiring branch leaders 
###################################


row1 <- round(c(mean(d$cxn_down_percentage_correct_full, na.rm=T), mean(d$cxn_down_percentage_correct_full, na.rm=T)*24), 2)
row2 <- round(c(mean(d$cxn_down_percentage_correct_full[d$urban==1], na.rm=T), mean(d$cxn_down_percentage_correct_full[d$urban==1], na.rm=T)*24), 2)
row3 <- round(c(mean(d$cxn_down_percentage_correct_full[d$urban==0], na.rm=T), mean(d$cxn_down_percentage_correct_full[d$urban==0], na.rm=T)*24), 2)
row4 <- row3-row2
row5 <- round(c(mean(d$faces_nonakan_percentage_correct_full, na.rm=T), mean(d$faces_nonakan_percentage_correct_full, na.rm=T)*24), 2)
row6 <- round(c(mean(d$faces_akan_percentage_correct_full, na.rm=T), mean(d$faces_akan_percentage_correct_full, na.rm=T)*24), 2)
row7 <- row6-row5
row8 <- round(c(mean(d$cxn_down_percentage_correct_full[d$old_branch_position==0], na.rm=T), mean(d$cxn_down_percentage_correct_full[d$old_branch_position==0], na.rm=T)*24), 2)
row9 <- round(c(mean(d$cxn_down_percentage_correct_full[d$old_branch_position==1], na.rm=T), mean(d$cxn_down_percentage_correct_full[d$old_branch_position==1], na.rm=T)*24), 2)
row10 <- row9-row8

names <- c("Full sample", "Urban", "Rural", "Difference", "Non Akans", "Akans", "Difference",  "Non-incumbents", "Incumbents", "Difference")

table_connections <- rbind(row1, row2, row3, row4, row5, row6, row7, row8, row9, row10)  
table_connections <- cbind(names, table_connections)
xtable(table_connections) 

### Relevant t-tests

ttesta <- t.test(d$cxn_down_percentage_correct_full[d$urban==1], d$cxn_down_percentage_correct_full[d$urban==0])

ttestb <- t.test(d$faces_nonakan_percentage_correct_full, d$faces_akan_percentage_correct_full)

ttestc <- t.test(d$cxn_down_percentage_correct_full[d$old_branch_position==1], d$cxn_down_percentage_correct_full[d$old_branch_position==0])


ttesta
ttestb
ttestc


###################################
###### Figure 6
###################################

#Over a third of branch leaders (35%) engaged in six or more of the nine campaign activities, while 28% performed three or less. 
summary(d$cpgn_brok_index[d$old_branch_position==1])
tab <- table(d$cpgn_brok_index[d$old_branch_position==1])
prop.table(tab)*100

#There is also significant variation across incumbent branch leaders in their post-campaign brokerage. While more than half of respondents engaged in none of the three post-election activities, more than one- third engaged in all three. 
summary(d$broker_up_down_index[d$old_branch_position==1])
tab2 <- table(d$broker_up_down_index[d$old_branch_position==1])
prop.table(tab2)*100

#CODE FOR FIGURE 6

reg1 <- lm(cpgn_brok_index ~ cxn_up_percentage_correct_full + cxn_down_percentage_correct_full  + aspires_higherposition + female  + log(age)+ local_eth_minority  + as.factor(bio_educ_three) + asset_index + years_comm + chief_relative    + old_ps_org  + old_ps_sec + old_ps_womenorg  + old_ps_youthorg+ ethfrac_upperlevel_2km  +branch_is_active.y + wealth_index_2km  + npp2016_pres_p.NEW + as.factor(constituency), data=d[d$old_branch_position==1,])


#POST-CAMPAIGN INDEX 
reg3 <- lm(broker_up_down_index ~ cxn_up_percentage_correct_full + cxn_down_percentage_correct_full  + aspires_higherposition + female  + log(age)+ local_eth_minority  + as.factor(bio_educ_three) + asset_index + years_comm + chief_relative    + old_ps_org  + old_ps_sec + old_ps_womenorg  + old_ps_youthorg+ ethfrac_upperlevel_2km  +branch_is_active.y + wealth_index_2km  + npp2016_pres_p.NEW + as.factor(constituency), data=d[d$old_branch_position==1,])


#### LEFT HAND SIDE PLOT OF FIGURE 6 

point_up <- reg1$coefficients[2]
segment_up1 <- reg1$coefficients[2]+1.96*summary(reg1)$coef[2,2]
segment_down1 <- reg1$coefficients[2]-1.96*summary(reg1)$coef[2,2]

point_up2 <- reg3$coefficients[2]
segment_up2 <- reg3$coefficients[2]+1.96*summary(reg3)$coef[2,2]
segment_down2 <- reg3$coefficients[2]-1.96*summary(reg3)$coef[2,2]

point_down <- reg1$coefficients[3]
segment_up3 <- reg1$coefficients[3]+1.96*summary(reg1)$coef[3,2]
segment_down3 <- reg1$coefficients[3]-1.96*summary(reg1)$coef[3,2]

point_down2 <- reg3$coefficients[3]
segment_up4 <- reg3$coefficients[3]+1.96*summary(reg3)$coef[3,2]
segment_down4 <- reg3$coefficients[3]-1.96*summary(reg3)$coef[3,2]

### PLOT 
par(mfrow=c(1,2))
par(mai=c(1.4,2.0,0.42,0.5))

points_x <- c(point_up2, point_up)
points_y <- c(1:2)
plot(points_x, points_y, xlim=c(-1,4), axes=F, pch=19, ylab="", xlab="Coefficient", col="dodgerblue", ylim=c(0, 3), cex=1.6)
abline(v=0, lty=2) 
points(c(point_down2, point_down), c(1.3, 2.3), pch=17, col="firebrick1", cex=1.6)
axis(1)
axis(2, labels=c(expression(bold("Post-campaign index")), expression(bold("Campaign index"))), at=c(1.1, 2.2), las=2)
legend(#x=1.4, y=0.1, 
  x=0.5, y=1,
  c("Connections down", "Connections up"), pch=c(17,19), col=c("firebrick1", "dodgerblue"), bty="n", cex=0.9) 
segments(segment_up2, 1, segment_down2, 1)
segments(segment_up4, 1.3, segment_down4, 1.3)
segments(segment_up3, 2.3, segment_down3, 2.3)
segments(segment_up1, 2, segment_down1, 2)

#### RIGHT HAND SIDE PLOT OF FIGURE 6 

## Calculate CNX_UP value at 1sd above and below mean
sd_above <- mean(d$cxn_up_percentage_correct_full[d$old_branch_position==1])+ sd(d$cxn_up_percentage_correct_full[d$old_branch_position==1])
sd_below <- mean(d$cxn_up_percentage_correct_full[d$old_branch_position==1])- sd(d$cxn_up_percentage_correct_full[d$old_branch_position==1])
# 
# ## Calculate CNX_DOWN value at 1sd above and below mean
sd_above2 <- mean(d$cxn_down_percentage_correct_full[d$old_branch_position==1])+ sd(d$cxn_down_percentage_correct_full[d$old_branch_position==1])
sd_below2 <- mean(d$cxn_down_percentage_correct_full[d$old_branch_position==1])- sd(d$cxn_down_percentage_correct_full[d$old_branch_position==1])
# 
## Create two new datasets -- one where cnx_up is 1sd above, and 1sd below mean
data_above <- data_below <- d[d$old_branch_position==1,] 
data_above$cxn_up_percentage_correct_full <- sd_above
data_below$cxn_up_percentage_correct_full <- sd_below
# 
# ## Create two new datasets -- one where cnx_down is 1sd above, and 1sd below mean
data_above2 <- data_below2 <- d[d$old_branch_position==1,] 
data_above2$cxn_down_percentage_correct_full <- sd_above2
data_below2$cxn_down_percentage_correct_full <- sd_below2
# 

## Store quantities of interest
res_mat <- matrix(NA,nrow = 12,ncol = 3)
## Loop over the 12 DVs 
dv <- c("helps_party_idcomm",
        "helps_comm_w_party",
        "helps_comm_w_da",
        "campaign_assistance",
        "campaign_jobs",
        "campaign_workwchief",
        "campaign_drivevoters",
        "campaign_orgtransport",
        "campaign_orgevents",  
        "campaign_orgrallies",
        "campaign_gifts",
        "campaign_canvass")
## Loop 
for (i in 1:12){
  fmreg <- as.formula(paste(dv[i],"~","cxn_up_percentage_correct_full + cxn_down_percentage_correct_full  + aspires_higherposition + female  + log(age)+ local_eth_minority  + as.factor(bio_educ_three) + asset_index + years_comm + chief_relative    + old_ps_org  + old_ps_sec + old_ps_womenorg  + old_ps_youthorg +ethfrac_upperlevel_2km + branch_is_active.y + wealth_index_2km  + npp2016_pres_p.NEW + as.factor(constituency)",sep = " "))
  ## Formula
  reg_a <- glm(fmreg,  data=d[d$old_branch_position==1,], family="binomial")
  ## Calculate predicted values using marg (in "modmarg" package)
  margins1 <- marg(reg_a, var_interest = 'cxn_up_percentage_correct_full', type = 'levels', at_var_interest = c(sd_above, sd_below))
  upper1 <- margins1[[1]][1,2]
  lower1 <- margins1[[1]][2,2]
  se_up  <- margins1[[1]][1,3]
  se_down <- margins1[[1]][2,3]
  ## Want the difference in the predicted probs
  diff <- upper1 - lower1 
  ## Need the se of the difference 
  se_diff <- sqrt(se_up^2 + se_down^2)
  ci_up <- diff+(1.96*se_diff)
  ci_down <- diff-(1.96*se_diff) 
  ## Things you want to plot/keep are: diff, ci_up, and ci_down
  res_mat[i,] <-c(diff, ci_up, ci_down)
}

############################
### Do the same for connections down 
############################

## Store quantities of interest
res_mat2 <- matrix(NA,nrow = 12,ncol = 3)
## Loop 
for (i in 1:12){
  ## Vars for regression 
  fmreg <- as.formula(paste(dv[i],"~","cxn_up_percentage_correct_full + cxn_down_percentage_correct_full  + aspires_higherposition + female  + log(age)+ local_eth_minority  + as.factor(bio_educ_three) + asset_index + years_comm + chief_relative    + old_ps_org  + old_ps_sec + old_ps_womenorg  + old_ps_youthorg +ethfrac_upperlevel_2km + branch_is_active.y + wealth_index_2km  + npp2016_pres_p.NEW + as.factor(constituency)",sep = " "))
  ## Formula
  reg_a <- glm(fmreg,  data=d[d$old_branch_position==1,], family="binomial")
  ## Calculate predicted values using marg (in "modmarg" package)
  ## Simulating 1 sd above and below the mean of the var of interest
  margins1 <- marg(reg_a, var_interest = 'cxn_down_percentage_correct_full', type = 'levels', at_var_interest = c(sd_above2, sd_below2))
  upper1 <- margins1[[1]][1,2]
  lower1 <- margins1[[1]][2,2]
  se_up  <- margins1[[1]][1,3]
  se_down <- margins1[[1]][2,3]
  ## Want the difference in the predicted probs
  diff <- upper1 - lower1 
  ## Need the se of the difference 
  se_diff <- sqrt(se_up^2 + se_down^2)
  ci_up <- diff+(1.96*se_diff)
  ci_down <- diff-(1.96*se_diff) 
  ## Things you want to plot/keep are: diff, ci_up, and ci_down
  res_mat2[i,] <-c(diff, ci_up, ci_down)
}

### Plotting the output

par(mai=c(1.4,3,0.42,1.0))

y_var <- c(1:15)
x_var <- seq(0, 0.20, by=0.014)

plot(x=x_var, y=x_var, xlim=c(-0.15, 0.20), pch=19, axes=F, ylab="", xlab="Difference in predicted probability", col="white", ylim=c(0, 18), cex=1.6, cex.axis=1.4, cex.lab=1.4)
dvs <- c(
  "Help party identify recipients",
  "Help voters contact party",
  "Help voters contact local govt.",
  expression(bold("Post-campaign index")), 
  "", 
  "Provide campaign assistance",
  "Find jobs for residents",
  "Work with the chief",
  "Drive voters to vote",
  "Organize transport for voters",
  "Organize community events",  
  "Organize rallies",
  "Distribute gifts",
  "House-to-house canvass",
  expression(bold("Campaign index"))
)
axis(2, labels=dvs, at=seq(from=1.1, to=15.1, by=1), las=2, cex.lab=1.5)
y_list <- c(1,2,3,6:14)
points(x=res_mat[,1], y=y_list, pch=19, col="dodgerblue")
abline(v=0, lty=2) ## Plot Connections up
segments(res_mat[,2], y_list, res_mat[,3], y_list, col="black", lty=1)
axis(1)

y_list2 <- c(1.2,2.2,3.2,6.2:14.2)
points(res_mat2[,1], y_list2, pch=17, col="firebrick1")
segments(res_mat2[,2], y_list2, res_mat2[,3], y_list2, col="black", lty=1)

legend(x=-0.015, y=0.6, c("Connections down", "Connections up"), pch=c(17,19), col=c("firebrick1", "dodgerblue"), bty="n", cex=0.9) 



###################################
###### Table 3: Mode of branch leader selection
###################################


table(ps$contestation_type)
table(ps$contestation_type) / sum(table(ps$contestation_type))


###################################
######  Table 4: Predictors of selection mode 
###################################

#Table 4 code
m1 <- glm(all_acclaimed ~ cxn_down_percentage_correct_full  + cxn_up_percentage_correct_full  + closest_ps_km  + wealth_index_2km  + npp2016_pres_p.NEW  + ethfrac_upperlevel_2km  + branch_is_active  + as.factor(constituency), data=bl, family="binomial")
summary(m1)

m2 <- glm(all_acclaimed ~ cxn_down_percentage_correct_full  + cxn_up_percentage_correct_full  + closest_ps_km  + wealth_index_2km  + npp12to16_ps_swing.NEW  + ethfrac_upperlevel_2km  + branch_is_active  + as.factor(constituency), data=bl, family="binomial")
summary(m2)

m3 <- glm(multiple_aspirants ~ cxn_down_percentage_correct_full  + cxn_up_percentage_correct_full  + incumbent_is_contesting  + branch_is_active.1 + wealth_index_2km + ethfrac_upperlevel_2km + npp2016_pres_p.NEW + female  + log(age)  + as.factor(bio_educ_three)  + asset_index  + as.factor(position.fac)  + as.factor(constituency), data= ps, family="binomial")
summary(m3)

m4 <- glm(multiple_aspirants~ cxn_down_percentage_correct_full + cxn_up_politicians  +  cxn_up_bureaucrats  +  cxn_up_constexecs + incumbent_is_contesting  + branch_is_active.1 + wealth_index_2km  + ethfrac_upperlevel_2km + npp2016_pres_p.NEW + female  + log(age) + as.factor(bio_educ_three)  + asset_index  + as.factor(position.fac)  + as.factor(constituency),  data= ps, family="binomial")
summary(m4)

m5 <- glm(backroom_deal ~ cxn_down_percentage_correct_full  +  cxn_up_politicians  +  cxn_up_bureaucrats  +  cxn_up_constexecs +  incumbent_is_contesting, data= ps[ps$multiple_aspirants==1,], family="binomial")
summary(m5)

m6 <- glm(backroom_deal ~ cxn_down_percentage_correct_full +  cxn_up_politicians  +  cxn_up_bureaucrats  +  cxn_up_constexecs +  incumbent_is_contesting + branch_is_active.1+ wealth_index_2km + ethfrac_upperlevel_2km + npp2016_pres_p.NEW + as.factor(position.fac), data= ps[ps$multiple_aspirants==1,], family="binomial")
summary(m6)

m7 <- glm(backroom_deal~ cxn_down_percentage_correct_full  + cxn_up_politicians  +  cxn_up_bureaucrats  +  cxn_up_constexecs  + branch_is_active.y + wealth_index_2km  + ethfrac_upperlevel_2km + npp2016_pres_p.NEW + as.factor(ps2018_contest_position), data= d.incrun[d.incrun $multiple_aspirants==1,], family="binomial")
summary(m7)



#Simulating from column 1, branches with incumbent branch leaders that have average connections down at the 90th percentile (49%) are 25.6 percentage points less likely (95% CI: 4.1, 43.9) to have all positions go unopposed than branches at the 10th per- centile (1%). 

m1 <- glm(all_acclaimed ~ cxn_down_percentage_correct_full  + cxn_up_percentage_correct_full  + closest_ps_km  + wealth_index_2km  + npp2016_pres_p.NEW  + ethfrac_upperlevel_2km  + branch_is_active  + as.factor(constituency), data=bl, family="binomial")
summary(m1)

#MAKE DUMMIES: 
bl$assin_south <- ifelse(bl $constituency=="assin_south", 1, 0)
bl $atwima_mponua <- ifelse(bl $constituency=="atwima_mponua", 1, 0)
bl $ayawaso_central <- ifelse(bl $constituency=="ayawaso_central", 1, 0)
bl $bortianor <- ifelse(bl $constituency=="bortianor", 1, 0)
bl $dome_kwabenya <- ifelse(bl $constituency=="dome_kwabenya", 1, 0)
bl $ho_central <- ifelse(bl $constituency=="ho_central", 1, 0)
bl $lower_manya_krobo <- ifelse(bl $constituency=="lower_manya_krobo", 1, 0)
bl $manhyia_south <- ifelse(bl $constituency=="manhyia_south", 1, 0)
bl $north_tongu <- ifelse(bl $constituency=="north_tongu", 1, 0)

quantile(bl$cxn_down_percentage_correct_full, c(.1, .9), na.rm=T)
      # 10%       90% 
# 0.0100000 0.4866667 

bl.dat <- na.omit(bl[, c("all_acclaimed" ,"cxn_down_percentage_correct_full" , "cxn_up_percentage_correct_full" , "closest_ps_km" , "wealth_index_2km" , "npp2016_pres_p.NEW" , "ethfrac_upperlevel_2km" , "branch_is_active" ,  "assin_south", "atwima_mponua", "ayawaso_central", "bortianor", "dome_kwabenya", "ho_central", "lower_manya_krobo", "manhyia_south", "north_tongu")])

set.seed(48103)
vcovCL <- vcov(m1)
N <- 1000
betas <- mvrnorm(n = N, mu= m1 $coefficients, vcovCL)

xt <- cbind(1, 0.4866667, bl.dat$cxn_up_percentage_correct_full , bl.dat$closest_ps_km, bl.dat$wealth_index_2km, bl.dat$npp2016_pres_p.NEW, bl.dat$ethfrac_upperlevel_2km, bl.dat$branch_is_active, bl.dat$assin_south, bl.dat$atwima_mponua, bl.dat$ayawaso_central, bl.dat$bortianor, bl.dat$dome_kwabenya, bl.dat$ho_central, bl.dat$lower_manya_krobo, bl.dat$manhyia_south, bl.dat$north_tongu)

ct <- cbind(1, 0.0100000, bl.dat$cxn_up_percentage_correct_full , bl.dat$closest_ps_km, bl.dat$wealth_index_2km, bl.dat$npp2016_pres_p.NEW, bl.dat$ethfrac_upperlevel_2km, bl.dat$branch_is_active, bl.dat$assin_south, bl.dat$atwima_mponua, bl.dat$ayawaso_central, bl.dat$bortianor, bl.dat$dome_kwabenya, bl.dat$ho_central, bl.dat$lower_manya_krobo, bl.dat$manhyia_south, bl.dat$north_tongu)

tprob2<-apply((1/(1+exp(-(as.matrix(xt))%*% t(betas)))),1,as.vector)
cprob2<-apply((1/(1+exp(-(as.matrix(ct))%*% t(betas)))),1,as.vector)

tprob2.means <- apply(tprob2, 1, mean)
cprob2.means <- apply(cprob2, 1, mean)

diffprob2 <- tprob2 - cprob2
diffprob2.means <- apply(diffprob2, 1, mean)
mean(diffprob2.means)
#-0.2559649
quantile(diffprob2.means, c(0.025, 0.975))
       # 2.5%       97.5% 
# -0.43933848 -0.04136029 


#Simulating from column 6, a “backroom deal” is 25.8 percentage points (95% CI: 1.3, 46.4) more likely for positions in which the incumbent branch leader has upward ties to constituency leaders at the 90th percentile (60%) than the 10th percentile (0%). 

m6 <- glm(backroom_deal ~ cxn_down_percentage_correct_full 
+  cxn_up_politicians  +  cxn_up_bureaucrats  +  cxn_up_constexecs 
+  incumbent_is_contesting 
+ branch_is_active.1
+ wealth_index_2km 
+ ethfrac_upperlevel_2km
+ npp2016_pres_p.NEW + as.factor(position.fac), data= ps[ps$multiple_aspirants==1,], family="binomial")
summary(m6)

#MAKE DUMMIES: 
ps$org_dummy <- ifelse(ps$position.fac=="org", 1, 0)
ps$sec_dummy <- ifelse(ps$position.fac=="secretary", 1, 0)
ps$womenorg_dummy <- ifelse(ps$position.fac=="women_org", 1, 0)
ps$youthorg_dummy <- ifelse(ps$position.fac=="youth_org", 1, 0)

quantile(ps $cxn_up_constexecs[ps$multiple_aspirants==1], c(.1, .9), na.rm=T)
# 10% 90% 
# 0.0 0.6 

ps.dat <- na.omit(ps[, c("backroom_deal" ,"cxn_down_percentage_correct_full" , "cxn_up_politicians" , "cxn_up_bureaucrats", "cxn_up_constexecs", "incumbent_is_contesting" ,  "branch_is_active.1",  "wealth_index_2km" , "ethfrac_upperlevel_2km",  "npp2016_pres_p.NEW", "org_dummy", "sec_dummy",  "womenorg_dummy", "youthorg_dummy")])

set.seed(48109)
vcovCL <- vcov(m6)
N <- 5000
betas <- mvrnorm(n = N, mu= m6 $coefficients, vcovCL)

xt <- cbind(1, ps.dat$cxn_down_percentage_correct_full , ps.dat$cxn_up_politicians , ps.dat$cxn_up_bureaucrats , 0.6, ps.dat$incumbent_is_contesting , ps.dat$branch_is_active.1 , ps.dat$wealth_index_2km , ps.dat$ethfrac_upperlevel_2km , ps.dat$npp2016_pres_p.NEW , ps.dat$org_dummy , ps.dat$sec_dummy , ps.dat$womenorg_dummy , ps.dat$youthorg_dummy)

ct <- cbind(1, ps.dat$cxn_down_percentage_correct_full , ps.dat$cxn_up_politicians , ps.dat$cxn_up_bureaucrats , 0, ps.dat$incumbent_is_contesting , ps.dat$branch_is_active.1 , ps.dat$wealth_index_2km , ps.dat$ethfrac_upperlevel_2km , ps.dat$npp2016_pres_p.NEW , ps.dat$org_dummy , ps.dat$sec_dummy , ps.dat$womenorg_dummy , ps.dat$youthorg_dummy)

tprob2<-apply((1/(1+exp(-(as.matrix(xt))%*% t(betas)))),1,as.vector)
cprob2<-apply((1/(1+exp(-(as.matrix(ct))%*% t(betas)))),1,as.vector)

tprob2.means <- apply(tprob2, 1, mean)
cprob2.means <- apply(cprob2, 1, mean)

diffprob2 <- tprob2 - cprob2
diffprob2.means <- apply(diffprob2, 1, mean)
mean(diffprob2.means)
#0.2579398
quantile(diffprob2.means, c(0.025, 0.975))
  #       2.5%      97.5% 
#0.01277301 0.46445201 


#Moreover, there is no evidence that the pool of aspirants is somehow pre-screened to weed out brokers with especially low downward ties. Indeed, the median branch leader protected through a backroom deal had very low connections down: 12.5%. This is lower than the median among all respondents (17%).

summary(d$cxn_down_percentage_correct_full [d$won_seat==1 & d$backroom_deal==1 &is.na(d$backroom_deal)==F & is.na(d$won_seat)==F])
summary(d$cxn_down_percentage_correct_full )




